The genetic structure and demographic history revealed by whole-genome resequencing provide insights into conservation of critically endangered Artocarpus nanchuanensis

Introduction Whole-genome resequencing technology covers almost all nucleotide variations in the genome, which makes it possible to carry out conservation genomics research on endangered species at the whole-genome level. Methods In this study, based on the whole-genome resequencing data of 101 critically endangered Artocarpus nanchuanensis individuals, we evaluated the genetic diversity and population structure, inferred the demographic history and genetic load, predicted the potential distributions in the past, present and future, and classified conservation units to propose targeted suggestions for the conservation of this critically endangered species. Results Whole-genome resequencing for A. nanchuanensis generated approximately 2 Tb of data. Based on abundant mutation sites (25,312,571 single nucleotide polymorphisms sites), we revealed that the average genetic diversity (nucleotide diversity, π) of different populations of A. nanchuanensis was relatively low compared with other trees that have been studied. And we also revealed that the NHZ and QJT populations harboured unique genetic backgrounds and were significantly separated from the other five populations. In addition, positive genetic selective signals, significantly enriched in biological processes related to terpene synthesis, were identified in the NHZ population. The analysis of demographic history of A. nanchuanensis revealed the existence of three genetic bottleneck events. Moreover, abundant genetic loads (48.56% protein-coding genes) were identified in Artocarpus nanchuanensis, especially in genes related to early development and immune function of plants. The predication analysis of suitable habitat areas indicated that the past suitable habitat areas shifted from the north to the south due to global temperature decline. However, in the future, the actual distribution area of A. nanchuanensis will still maintain high suitability. Discussion Based on total analyses, we divided the populations of A. nanchuanensis into four conservation units and proposed a number of practical management suggestions for each conservation unit. Overall, our study provides meaningful guidance for the protection of A. nanchuanensis and important insight into conservation genomics research.

Introduction: Whole-genome resequencing technology covers almost all nucleotide variations in the genome, which makes it possible to carry out conservation genomics research on endangered species at the whole-genome level.
Methods: In this study, based on the whole-genome resequencing data of 101 critically endangered Artocarpus nanchuanensis individuals, we evaluated the genetic diversity and population structure, inferred the demographic history and genetic load, predicted the potential distributions in the past, present and future, and classified conservation units to propose targeted suggestions for the conservation of this critically endangered species.
Results: Whole-genome resequencing for A. nanchuanensis generated approximately 2 Tb of data. Based on abundant mutation sites (25,312,571 single nucleotide polymorphisms sites), we revealed that the average genetic diversity (nucleotide diversity, p) of different populations of A. nanchuanensis was relatively low compared with other trees that have been studied. And we also revealed that the NHZ and QJT populations harboured unique genetic backgrounds and were significantly separated from the other five populations. In addition, positive genetic selective signals, significantly enriched in biological processes related to terpene synthesis, were identified in the NHZ population. The analysis of demographic history of A. nanchuanensis revealed the existence of three genetic bottleneck events. Moreover, abundant genetic loads (48.56% protein-coding genes) were identified in Artocarpus nanchuanensis, especially in genes related to early development and immune function of plants. The

Introduction
The Dalou Mountains, with the famous Jinfo Mountains as the highest peak, are the boundary between the Guizhou Plateau and the Sichuan Basin (Loṕez-Pujol et al., 2011;Deng, 2019;Zu et al., 2022;Zu et al., 2023;Supplementary Figure 1). Local orogeny led to the formation of its present complex topography, which is mainly composed of low mountain canyons and middle mountain platforms. The humid and warm climate and high habitat heterogeneity together contributed to the abundant biodiversity of the Dalou Mountains. A recent study indicated that there were more than 5,000 species of vascular plants in the Dalou Mountains, including more than 200 plants for which type specimens were collected from the Dalou Mountains (Ding et al., 2014;Deng, 2019). Moreover, due to the Qinling Mountains-Daba Mountains in the north, the Dalou Mountains were less affected by the Quaternary ice age and retained many ancient relict species, such as Ginkgo biloba, Cathaya argyrophylla and Taxus chinensis (Wang and Ge, 2006;Tang et al., 2012;Qian et al., 2016;Supplementary Figure 1). Although the Dalou Mountains harbour abundant biodiversity, it has been given less attention than Daba Mountains and Wuling Mountain, which are also located in Chongqing City and designated as priority areas for biodiversity conservation. This has greatly hindered the further protection of endangered plants in this region, especially species with small populations.
Currently, many endangered and narrowly distributed plant species are found in the Dalou Mountains (Deng, 2019). For the rescue of endangered plants, ex -situ protection is generally considered to be an important and reliable means to prevent extinction (Huang and Zhang, 2012;Xu and Zang, 2022). However, recent studies have shown that ex -situ conservation of endangered plants with a lack of comprehensive population genetic information often leads to failure to cover the genetic diversity of wild populations and achieve protection of germplasm resources (Wu et al., 2013;Xu and Zang, 2022). Therefore, it is very important to identify conservation units for the protection and management of endangered species, which could provide a clear understanding of the distribution of genetic diversity and the boundary of population units for protection managers (Funk et al., 2012;Li et al., 2020;. With recent developments in sequencing technology and reductions in sequencing costs, whole-genome resequencing has become increasingly popular because it covers almost all nucleotide variations in the genome (Hohenlohe et al., 2021;Ma et al., 2021;Ma et al., 2022). This has greatly helped accurately and clearly reveal the demographic history and genetic structure of endangered species, which are further helpful for the division of conservation units and the development of effective conservation strategies (Ma et al., 2021;Ma et al., 2022;. Climate change is considered to be the major threat to biodiversity in the 21 st century because it may lead to loss of biodiversity, termination of evolutionary potential and disruption of ecological services (Thomas et al., 2004;Pimm, 2008;Dawson et al., 2011). For endangered species that are vulnerable to climate fluctuations, it is difficult to formulate appropriate conservation plans without considering the impact of climate change (Huang, 2011). Recent studies have indicated that climate change greatly threatens the population maintenance of endangered species and also drive the distribution pattern of endangered species to change, affecting the effectiveness of species conservation (Ackerly et al., 2010;Barber et al., 2016;Xue et al., 2022). Species distribution modelling combines species distribution records with environmental variables to predict the potential distribution range of species and has been widely used in the assessment of biodiversity conservation and priority management of many endangered species (Qin et al., 2017a;Tang et al., 2018;Gao et al., 2022;Xue et al., 2022).
Artocarpus nanchuanensis S.S. Chang, S.C. Tan, Z.Y. Liu, which is mainly distributed in the northern section of the Dalou Mountains, is the northernmost species in the natural distribution of the genus Artocarpus (Wu and Zhang, 1989; Figure 1; Supplementary Figure 1). This species is distinguished from other species of Artocarpus by two rows of alternate leaves, elliptic to nearly round leaves, irregular spherical polythalamic fruit with papillary protrusion on the surface, etc (Figure 1). Recent phylogenetic studies have indicated that this species belongs to the basal clade of the subg. Pseudojaca (Ser. Clavati), which contains mostly Chinese species (Gardner et al., 2021). The fruit of A. nanchuanensis is rich in nutrients and can be used for brewing, processing jam and beverages (He, 2014;He et al., 2022a). Straight trunks are often used for shipbuilding and making furniture (He, 2014). Although artificial breeding technologies for A. nanchuanensis have gradually been established in recent years owing to its high potential application value, only breeding bases and several parks currently have cultivation (Tan et al., 2015). In addition, due to its fragmented distribution and small wild population, this species was rated as critically endangered (CR) and as a Wild Plants with Extremely Small Populations (WPESP) in Chongqing City (Qin et al., 2017b). A recent field investigation indicated that A. nanchuanensis was distributed only in 4 counties of Chongqing City, and no more than 100 maternal trees remained (He, 2014). It was also found that adult trees were rare in most populations, and there was a pattern of more seedlings and fewer young trees (field survey data). Therefore, research on conservation biology and molecular genetics for this critically endangered species is urgently needed to support the formulation of further conservation plans.
The recently reported genome assembly at the chromosome level of A. nanchuanensis provided an opportunity to resequence at the whole genome level (He et al., 2022a). To reveal the genetic structure and demographic history of A. nanchuaensis fully and accurately and provide meaningful guidance for its conservation, 101 individuals from 7 populations of A. nanchuanensis were resequenced in this study. We aimed to: (1) evaluate the genetic diversity and population structure of A. nanchuanensis; (2) infer the demographic history of this species and evaluate its genetic load; (3) predict the potential distributions of this species in the past, present and future; and (4) use our analysis results to classify conservation units and propose targeted suggestions for the conservation of this critically endangered species.
2 Materials and methods 2.1 Sample collection and whole genome resequencing   and then sampled seedlings at intervals of three meters. The number of individuals sampled from each population ranged from 11 to 27. Two populations from Yongchuan District were merged due to following reasons. (1) These two populations are located in the same branch range of Dalou Mountains (Supplementary Figure 1).
(2) Poorly differentiated and similar genetic backgrounds were detected between them.
(3) These two are closest, but far away from populations in other regions. A tender leaf of each selected individual was collected and taken back to the laboratory for storage at -80 degrees Celsius. Genomic DNA was extracted from samples using the standard cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1987). A Qubit4 Fluorometer nucleic acid protein fluorescence quantitative instrument and agarose gel electrophoresis were used to detect the concentration and quality of total genomic DNA, respectively. To construct paired-end libraries, we performed following treatments on the total genomic DNA (1 mg) of each sample: fragmentation, end repairing, adaptor ligation, PCR amplification and PCR product cyclization . The Agilent 2100 Bioanalyzer (Agilent DNA 1000 Reagents) was used to detect the fragment size and concentration of the constructed library. All libraries were sequenced on the BGISEQ-500 sequencing platform, which uses optimized Combinatorial Probe-Anchor Synthesis (cPAS) and DNA nanospheres (DNB) core sequencing technology, producing 150 bp pair end reads .

Genome mapping and SNP calling
Raw data with adapter sequences or low-quality sequences were filtered using SOAPnuke with the following parameters: -n 0.01 -l 20 -q 0.3 -adaMR 0.25 -ada_trim -polyX 50 -minReadLen 150 . BWA-MEM (Li, 2013) and default parameters were used to map paired-end clean reads to the reference genome of A. nanchuanensis (He et al., 2022b). SAMtools (Danecek et al., 2021) and Picard (http:// picard.sourceforge.net/) were used to sort and delete duplicate reads, respectively. Then, we employed the mpileup + call (-multiallelic-caller -variants-only) and filter modules (-g 3 -G 10 -e "%QUAL<10") in BCFtools (Danecek et al., 2021) to call and filter variable sites in each sample, and the view module (-types snps) in BCFtools was used to extract single nucleotide polymorphisms sites (SNPs) to obtain Dataset 1, including 25,312,571 SNPs. To ensure high-quality SNPs, VCFtools v0.1.17 (Danecek et al., 2011) was utilized to further filter Dataset 1, with the following settings: -maf 0.01 -max-missing 0.9 -min-alleles 2 -max-alleles 2. After filtering, a total of 18,906,649 SNPs (Dataset 2) remained for downstream analysis. SnpEff 5.0e (Cingolani et al., 2012) was used for SNP annotation. We first used SnpEff 5.0e to construct a new database based on the assembly and annotation file of the A. nanchuanensis genome, and then annotated Dataset 1 by using the constructed database and default parameters.

Population genetic structure
Principal component analysis (PCA) was conducted by using Plink (Purcell et al., 2007), and the R package (ggplot2; Wickham, 2016) was used for visualization of the analysis results. Population structure was analyzed by Admixture (Alexander et al., 2009), which could infer the stratification relationships of populations. First, PLINK (Purcell et al., 2007) was used to filter the linkage sites and obtain a temporary file containing the sites that needed to be preserved from Dataset 2 with the following parameters: indep pairwise 50 10 0.2. Independent SNPs (Dataset 3) were extracted based on the temporary file and Dataset 2, yielding a total of 5,146,515 SNPs. Then, the bed file converted from the vcf file was used as the input file of Admixture for population structure analysis. According to the total number of natural populations, the number of clusters K was set from 2 to 7. The optimal K value was determined based on the actual distribution information and PCA results. Finally, the results were visualized by Poppelper (Francis, 2017) based on the Q matrix obtained from the Admixture analysis.

Phylogenetic analysis and estimation of population genetic parameters
Based on the latest research of Artocarpus, we downloaded the resequencing data of five other species of Artocarpus (A. hypargyreus, A. styracifolius, A. tonkinensis, A. petelottii, and A. gomezianus) to explore the phylogenetic status of A. nanchuensis (Gardner et al., 2021;Lin et al., 2022;Supplementary Table 2). The first four species are thought to be in the same clade as A. nanchuensis (ser. Clavati), and the last one belongs to the sisters clade (ser. Peltati) of ser. Clavati (Gardner et al., 2021). These data were processed with the method described in section 2.2 to obtain SNP data. This dataset, along with Dataset 3 of A. nanchuanensis, was used to construct the phylogenetic tree by IQ-TREE (Nguyen et al., 2015). We first converted the vcf file of the combined dataset to phylip format and then used the following settings for tree building: automatic model selection, unpartitioned, and 1000 fast bootstrap analyses. The obtained phylogenetic tree was rooted with A. gomezianus (Gardner et al., 2021). In addition, based on the annotation results from SnpEff 5.0e (Cingolani et al., 2012), we calculated the genetic diversity of synonymous (p S ) and nonsynonymous sites (p N ) of each population using VCFtools with a 100-kb sliding window and a step size of 10 kb. Genetic diversity was defined as the average number of pairwise nucleotide differences (nucleotide diversity, p) in this study (Nei and Li, 1979;Tajima, 1983;Corbett-Detig et al., 2015). VCFtools v0.1.17 (Danecek et al., 2011) was also used to calculate other population genetic parameters with a 100-kb sliding window and a step size of 10 kb based on Dataset 2, including fixation statistics (F ST ) and p ratio values between each population pair.

Selective sweep analysis
To identify the potential selective sweeps in the two special populations of A. nanchuanensis (NHZ and QJT), we selected the windows with high genetic differentiation (top 1% of F ST and top 1% of nucleotide diversity ratio values (p QJT /p NHZ )) as the candidate divergent regions (CDRs) by using Microsoft Excel. Candidate genes were obtained from CDRs and GO enrichment analysis was carried out by using TBtools . Then we extracted the main genes from the significantly enriched GO items and analyzed the clustering of these genes in the corresponding gene family. The complete members of the gene family in A. nanchuanensis were searched and confirmed by HMMER 3.0 (Finn et al., 2015) and the NCBI Conserved Domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). In addition, the genes of the corresponding gene family in Arabidopsis thaliana were downloaded from the Phytozome database (Goodstein et al., 2012) and used together with the genes of A. nanchuanensis to construct a phylogenetic tree.

Demographic history inference
To reconstruct the detailed demographic history of A. nanchuanensis, BCFtools (mpileup + call module) (Danecek et al., 2021) was used to obtain the consensus sequence between the sample and reference genome, and the fq2psmcfa module embedded in paired sequential Markov condensation (PSMC) (Li and Durbin, 2011) was used to convert its output file into fasta-like format. Then, the demographic history of each population was inferred by a PSMC model (Li and Durbin, 2011) with the following parameter settings: -N25 -t15 -r5 -p "4 + 25 * 2 + 4 + 6. The generation time was set to 11 years on the basis of years of field observation of A. nanchuanensis, and the mutation rate per generation was set to 2.86e-8 (11* 2.6e-9) (Gaut et al., 1996).

Characteristics of deleterious mutations
It is predicted that deleterious mutations (genetic load) will destroy gene function, which would significantly reduce the fitness of individuals in the population (Mattila et al., 2012). Moreover, loss of function (LOF) was the major type of deleterious mutations. To estimate the genetic load of A. nanchuanensis, we first defined the dominant alleles in all individuals (more than 50% of individuals had homozygous alleles) as the ancestral allele genotype (Hu et al., 2021). Then, the ancestral allele loci of each individual were extracted and annotated by using VCFtools v0.1.17 (Danecek et al., 2011) and SnpEff 5.0e (Cingolani et al., 2012), respectively. For each individual, we calculated the number of homozygous LOF sites and the number of genes containing one or more homozygous LOF sites divided by the total number of protein-coding genes (PCGs). Finally, gene ontology (GO) analysis was applied to characterize the detected genes containing homozygous LOF sites based on the genome annotation results of recent studies by using TBtools He et al., 2022a).

Ecological niche modelling
To further reveal the change in the distribution pattern of A. nanchuanensis, MaxEnt (version 3.3) (Phillips et al., 2006) was used to predict the past (mid Pliocene warm period -ca. 3.205 Ma, and MIS19-ca. 787 ka), current  and future (2070 years, RCP4.5) potential suitable habitat areas. Since five bioclimatic variables were missing in the two past periods, only 14 bioclimatic variables were downloaded for the prediction (Supplementary Table 3). The 14 bioclimatic variables of current climate data and paleoclimate data were downloaded from the Paleoclim database (Brown et al., 2018) with 2.5 min resolution. The future climate data (2070, BCC-CSM1-1 model, RCP 45) were downloaded from the WorldClim database (https://worldclim.org/ data/v1.4/cmip5_2.5m.html). To minimize the impact of multicollinearity and overfitting on the stability and quality of the model, the 14 bioclimatic variables were analyzed with the Pearson correlation test. When the correlation coefficient was > | 0.80|, one of the two variables was removed. Ultimately, six bioclimatic variables were selected to predict the suitable habitat areas for four time periods, including BIO1 (annual mean temperature), BIO4 (temperature seasonality), BIO8 (mean temperature of wettest quarter), BIO12 (annual precipitation), BIO14 (precipitation of driest month), and BIO15 (precipitation seasonality; Supplementary Table 3). We obtained 23 distribution records from the field investigations, which covered most of the known living areas of A. nanchuanensis (Supplementary Table 4). We also removed duplicate records in the same grid cell at a spatial resolution of 2.5 arc minutes to reduce spatial autocorrelation (Boria et al., 2014). Finally, a total of 18 distribution records were used in the analysis. In addition, we extracted the six bioclimate variables from these 18 records for principal component analysis (R package). The first two principal components (PC1 and PC2) were extracted from the analysis to estimate the main changes in bioclimatic variables represented by climate instability.

Whole-genome resequencing
The whole-genome resequencing of 101 individuals generated approximately 2 Tb of data. The average coverage, average depth and primary mapping rate of all samples were 98.00%, 39.23 and 97.57%, respectively (Supplementary Table 5). The annotation results of all SNPs indicated that 60.80% of the SNPs were located in intergenic regions, 13.12% were located in exon regions, and 8.18% were located in intron regions (Supplementary Table 6).

Population structure
The PCA indicated that the first three principal components (PCs) explained 16.2%, 12.9% and 10.7% of the total genetic variance, respectively. PC1, PC2, and PC3 significantly separated the NHZ and QJT populations from all other populations, but did not separate the remaining five populations (Figure 2A; Supplementary Figures 2, 3). The analysis of population structure showed that for the optimal K value (K=4), all samples were divided into four groups, including QJT, NHZ, QJS and four other populations ( Figure 2C). It is worth noting that the genetic backgrounds of the NHZ and QJT populations were quite simple, with only a few genetic mixings compared to other populations ( Figure 2C). The QJS population, which was geographically close to the QJT population (only separated by a 30-meter-wide river), had a distinct genetic mixture of QJT and the other four populations (NQQ, NSF, WS and YC). However, there was no significant difference in genetic background among the four populations from Nanchuan District, Wansheng and Yongchuan District (NQQ, NSF, WS and YC) ( Figure 2C). The F ST analysis indicated that the degrees of differentiation between various populations of A. nanchuanensis were relatively low, all less than 0.01 (Supplementary Table 7).

Phylogenetic and genetic parameter estimation
Phylogenetic analysis showed that A. nanchuanensis and four species of ser. Clavati (A. tonkinensis, A. petelotii, A. hypargyreus, and A. styracifolius) gathered together and formed a sister branch ( Figure 2B). Consistent with the PCA results, the NHZ and QJT populations were significantly separated from the other populations in the phylogenetic tree, while the remaining five populations (MIX) were mixed to varying degrees ( Figure 2B). In addition, the average genetic diversity of different populations of A. nanchuanensis was relatively low, with p S values ranging from 1.3×10 -3 to 1.33×10 -3 and p N values ranging from 2.26×10 -3 to 2.32×10 -3 , and there were no significant differences between populations (Supplementary Table 8).

Signals of selective sweeps between populations
In the selective sweeps analysis, we detected 376 genes of 375 CDRs that were under clear positive selection in the NHZ population (Supplementary Table 9; Supplementary Table 10).  Table 11). The phylogenetic tree constructed based on 57 detected TPS genes of A. nanchuanensis and 33 TPS genes of Arabidopsis thaliana indicated that these five TPS genes belong to the TPS-b subfamily, and they were notably clustered into a small clade with high bootstrap support (bootstrap percentage 100%; Figure 3C).

Demographic history of A. nanchuanensis
The results of PSMC analysis indicated that A. nanchuanensis experienced three genetic bottleneck events in its evolutionary history ( Figure 4A). The effective population size (N e ) of A. nanchuanensis began to decline slowly from 10 Ma and encountered the first genetic bottleneck at 2.82 Ma. From 2.2-0.2 Ma, the N e gradually expanded from 5.4x10 4 to 1.12x10 5 . However, the N e has been declining since 0.19 Ma, except for a small increase at 0.09 Ma. Except for the earlier first genetic bottleneck, this species experienced the second and third genetic bottlenecks at 0.14 Ma and 0.06 Ma, respectively. Finally, the N e of A. nanchuanensis was stable in the range of approximately 6.5 x10 4 -8.3x10 4 .

Estimation of deleterious mutations
Loss of function in a homozygous state could be used as an indicator of species fitness. The statistical analysis of the annotation results for each individual revealed no significant differences in genes containing homozygous LOF sites/PCGs between populations (Supplementary Table 12). The mean value for all individuals was 48.56%. A total of 285 GO terms were assigned in the enrichment analysis of genes containing homozygous LOF sites, including 220, 34 and 31 GO terms belonging to the biological process, molecular function and cell component categories, respectively (Supplementary Table 13). In the biological process category, multiple GO terms related to plant development in early stages, as well as potassium ion homeostasis (GO: 0055075) and tyrosine catalytic process (GO: 006572), were significantly enriched (corrected p value < 0.05, Supplementary Figure 6).

Niche modelling
The PCA of bioclimate variables showed that the first two PCs (PC1 and PC2) explained 88.5% of the total climate change variation experienced by A. nanchuanensis ( Figure 4B; Supplementary Figures 7, 8). Among the six bioclimatic variables, mean temperature of the wettest quarter (BIO8) and annual average temperature (BIO1) contributed most significantly to the first two PCs, followed by annual precipitation (BIO12) and precipitation seasonality (BIO15) ( Figure 4B; Supplementary Figure 9). The prediction results showed that the area under the operating characteristic curve (AUC) of the model used in the analysis was 0.998, indicating that the prediction model had reached an excellent threshold (Swets, 1988). This suggested that the model employed in this study could accurately describe the relationship between species and climate. The results of the prediction analysis showed that in the mid-Pliocene warm period (3.205 Ma), the suitable habitat areas of A. nanchuanensis were mainly distributed in central and southeastern China ( Figure 5A). However, the suitable habitat areas were predicted to shift from northern to southern China and were mainly distributed in the south parts of China in MIS19 (Pleistocene, ca. 787 ka; Figure 5B). The prediction results also indicated that the current suitable habitat areas, which were widely distributed in central and south-eastern China, were more extensive than the actual distribution ( Figure 5C). Although the future suitable habitat areas were predicted to shrink slightly compared with the current suitable habitat areas, the actual distribution area of A. nanchuanensis (Chongqing City) will still maintain high suitability ( Figures 5C, D).

Discussion
In this study, we sampled 101 individuals of A. nanchuanensis, covering most of the wild distribution areas, and revealed its complex population structure, demographic history, and genetic threats it currently faces by using whole-genome resequencing data. In addition, we identified several key factors that led to the formation of the distribution pattern of A. nanchuanensis by combining the analysis of niche models and demographic history. Our findings provide important guidance for meaningful conservation actions for this critically endangered WPESP.
Compared to that of other trees that have been studied, the genetic diversity of A. nanchuanensis was lower (Corbett-Detig et al., 2015;Hazzouri et al., 2015;Wang et al., 2016;Lin et al., 2022), which may be due to the following factors: (1) Previous studies have shown that genetic bottlenecks may lead to a large loss of genetic diversity due to enhanced genetic drift and a reduction in natural selection efficiency (Jangjoo et al., 2016;Ma et al., 2021;Ma et al., 2022). This seemed to be reasonable in the case of A. nanchuanensis because the demographic history showed three potential bottleneck events, which may have led to a significant reduction in its genetic diversity.
(2) The actual N e found in the field investigation (<200) was significantly smaller than the N e (>65,000) inferred from demographic history, indicating that the population of A. nanchuanensis has recently encountered strong pressure. During our field investigation in 2019, three small populations of A. nanchuanensis were found in Qianqiu village (Nanchuan district). When we visited the site again in 2022, we found that only a few young trees remained in one of the populations. The main reason was that its straight trunks of A. nanchuanensis are often used by local people for shipbuilding and furniture. Based on the inferred large N e value and the evidence that multiple existing populations in the four districts shared very similar genetic backgrounds, we speculated that A. nanchuanensis could be widely distributed in the mountains around the four districts. However, this may not be the case, possibly because human activities such as deforestation, road construction, and the expansion of farmland have drastically reduced the N e of A. nanchuanensis. Human interference and habitat fragmentation will further accelerate the impact of genetic drift, thereby reducing genetic diversity (Ma et al., 2021). (3) Studies on other Artocarpus species indicated that birds and large mammals may be potential seed disseminators (Dickinson et al., 2020). However, considering that most of the existing populations of A. nanchuanensis are located near settlement, the likelihood of the occurrence of large mammals is extremely low. Moreover, field investigations showed that mature individuals of A. nanchuanensis are rare (numbers between 1 and 5) in most populations, and the geographical distances between populations usually exceed 10 kilometers. These factors may together hinder the production of outbred seeds, which may be an important way to increase the genetic diversity of the population.
Artocarpus is a genus mainly distributed in tropical regions, and its diversity center is located in Borneo (Gardner et al., 2021). However, as the northernmost species of this genus, the reason for the formation of the distribution pattern of A. nanchuanensis has always been a mystery. Our study revealed that Quaternary climate change played an important role in forming the current distribution pattern of this species, especially temperature. A recent study indicated that the global temperature has been falling since 10 Ma (Westerhold et al., 2020). Combining the prediction analysis and demographic history, we speculated that due to the decrease of temperature in the early stage (10-2.82 Ma), the suitable habitat areas of A. nanchuanensis shrank, accompanied by a slow decline in N e (Figures 4A, 5A, B). With the advent of the Quaternary Ice Age, most areas of China were covered by glaciers, except for the mountainous areas in central and southern China, which harboured diverse terrain and a suitable climate, providing longterm stable habitats for plants (Axelrod et al., 1998;Loṕez-Pujol et al., 2011;Tang et al., 2012). Thus, A. nanchuanensis might have undergone migration and expansion in this region (south of the Qinling-Daba Mountains) during most of the Pleistocene (2.2-0.2 Ma) (Figures 4A, 5B). However, N e began to decline rapidly at 0.19 Ma, which coincided with the penultimate glaciation (0.19-0.14 Ma) (Colleoni et al., 2016). In general, during the Pleistocene glaciation (2000-12ka), the significant increase and decrease of Ne detected in A. nanchuanensis was consistent with the published research on other Spermatophytes and Ginkgo biloba, which is located in the same refuge -Dalou Mountains (Qiu et al., 2011;Zhao et al., 2019;Ma et al., 2021;Ma et al., 2022). In addition, the recent comparative analysis of demographic history of Artocarpus species showed that bioclimatic change and habitat were important factors in shaping the population history (Patil et al., 2022). Therefore, we inferred that similar habitat -acid soil of A. nanchuanensis and A. altilis may partly explain their highly consistent demographic histories (Patil et al., 2022).
The prediction results showed that the current suitable habitat areas of A. nanchuanensis were far larger than its actual distribution area ( Figure 5C), which may be due to the rapid climate fluctuations in the late Quaternary (approximately 0.8 -0.1 Ma, Hofreiter and Stewart, 2009), leading to the disappearance of populations in other areas, while only the populations in the refuge could survive. Previous studies have shown that the Dalou Mountains are an important refuge for plants, harbouring many living fossil plants, such as Ginkgo biloba, Cathaya argyrophylla, and Davidia involucrata. (Wang and Ge, 2006;Tang et al., 2012). Therefore, we speculated that only the remaining population of A. nanchuanensis located in Dalou Mountain survived in the face of the considerable challenges of climate change. In addition, little is known about the pollination and diffusion mechanisms of A. nanchuanensis, and further research is urgently needed in the future to comprehensively infer the underlying reasons for the differences between the current suitable habitat areas and actual distribution area. It is also worth noting that our prediction analysis did not include environmental variables, such as topographical variables and elevation, which should be integrated to yield more accurate prediction results in the future.
The GO annotation of positive selection signal regions of A. nanchuanensis indicated that the NHZ population was significantly enriched in 5 TPS genes, which were detected to have a close genetic relationship (Figure 3). TPS is the key enzyme responsible for terpene biosynthesis (Yu et al., 2020;Jia et al., 2022), which plays an important role in the defense mechanism of plants against herbivores and pathogenic microorganisms (Kessler, 2001;Gershenzon and Dudareva, 2007;Richter et al., 2015) and in the process of attracting pollinators (Yu et al., 2020;Zhao et al., 2021). There are only two adult trees and dozens of annual seedlings in the NHZ population, one of which is a listed ancient tree (DBH 63 cm), and the other adult tree (DBH 25 cm) is inferred as its offspring. This population is relatively isolated and now surrounded by cities (located within the crematorium area in Nanchuan District). Combined with the unique genetic background of the NHZ population revealed by the population structure analysis, it is inferred that the population has undergone a certain degree of differentiation during the long-term isolated growth process, especially in defense-related genes.
Our study also showed that there is a large amount of genetic load in A. nanchuanensis (Hu et al., 2020;Yi et al., 2022;Supplementary Table 12). Moreover, the annotation of genes containing homozygous LOF sites of all populations revealed that they were significantly enriched in multiple GO terms related to plant development in the early stage and plant immunity, such as potassium ion homeostasis (GO: 0055075) and tyrosine catabolic process (GO: 006572; Supplementary Figure S6). Previous studies have shown that potassium (K+) plays an important role in many basic processes in plants, including cell homeostasis, membrane transport, osmotic regulation and the immune response (Amtmann et al., 2008;Wang and Wu, 2015;Shi et al., 2018). Tyrosine phosphorylation has also been shown to be crucial for regulating immune signals in plants (Lin et al., 2014;Macho et al., 2014;Macho et al., 2015;Liu et al., 2018). All of these GO terms were related to deleterious mutations, which may reduce the ability of A. nanchuanensis to adapt to unsuitable environments. Field observations indicated that most populations included numerous newly sprouted seedlings, but perennial seedlings or saplings were rare. Therefore, it could be inferred that the seed germination rate of A. nanchuanensis is high, but the seedling survival rate is low, which may be related to the accumulation of deleterious mutations in its genes about development or immunity. However, further analysis of these deleterious mutations is needed to reach a definite conclusion.
Most of the populations of A. nanchuanensis are located outside nature reserves. In addition, the low genetic diversity and high genetic load, as well as serious human destruction and habitat disturbance, mean that the future survival of A. nanchuanensis may face great challenges. Based on our study results and actual conservation status, we propose to divide all populations into four conservation units, NHZ, QJT, QJS and four other populations (WS, NSF, NQQ and YC), and implement the following protection actions as soon as possible: (1) PCA, genetic structure and phylogenetic analysis revealed the unique genetic backgrounds of QJT and NHZ, which were different from those of other populations (Figures 2A, C). To protect the genetic diversity of A. nanchuanensis, priority should be given to the protection of these two small populations. The geographical distance between QJT and QJS was only 1 km, with the populations separated by a river approximately 30 meters wide. Surprisingly, our analysis showed that QJT and QJS had significant similarities and differences in genetic composition at the same time, which seemed to be an excellent example of geographical isolation and differentiation (Figures 2A, C). Therefore, we suggest establishing plant microreserves for these two small populations in Qijiang District, which harbour many apparent advantages, including lower costs, fast construction, and strong pertinence, and could also improve the public awareness of biodiversity conservation of A. nanchuanensis among local residents and even society (Fos et al., 2014;Laguna et al., 2016;Zhang, 2017). The NHZ population is restricted to flower beds of the crematorium in Nanchuan District, and one of the adult trees is listed as an ancient tree. Therefore, consulting with the management personnel of the crematorium is recommended to increase fencing for the flower beds, allowing this population to grow and breed naturally without human interference.
(2) Our field investigation revealed that the individual with the largest DBH (65 cm vs. 63 cm of the listed ancient tree from NHZ) was located in Honglu Town, Yongchuan District, and there is a sand factory under construction around the ancient individual. It is recommended to contact the local forestry administration to confirm the age of the ancient tree as soon as possible and apply in situ conservation for this small population with only three mature trees. (3) Considering the rarity of adult trees, isolated populations and small F ST values between various populations of A. nanchuanensis, artificial cross-pollination among four populations with similar genetic backgrounds should be carried out to maintain the genetic diversity and increase the survival potential of this species. (4) Based on the current status of the extremely small population of A. nanchuanensis, it is strongly recommended to pursue ex situ protection of the populations of the four conservation units to comprehensively preserve the germplasm resources and genetic diversity of this species. Considering the current and future habitat suitability of A. nanchuanensis, large-scale botanical gardens in central China are the preferred choices, such as the Nanshan Botanical Garden in Chongqing City and the Wuhan Botanical Garden of the Chinese Academy of Sciences in Hubei Province.

Conclusion
Based on the whole-genome resequencing data of 101 critically endangered A. nanchuanensis individuals, we revealed that repeated bottleneck events, human disturbance and habitat fragmentation, as well as excessive self-pollination within populations, together contributed to low genetic diversity. Quaternary climate change, especially temperature change, plays an important role in forming the current distribution pattern and endangered status of the species. The Dalou Mountains provided an important refuge for A. nanchuanensis during Quaternary climate fluctuations. Genetic analysis revealed that the NHZ population not only had a special genetic background but also had specific differentiation in defenserelated genes. The analysis also indicated that this species has an abundant genetic load, especially in the genes related to the early development and immune function of plants. Based on the study results, we divided the populations of A. nanchuanensis into four conservation units and proposed a number of practical management suggestions for each conservation unit of this critically endangered species.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI, PRJNA948357 (SRA).

Author contributions
HD designed the research. CX, BW and JZ collected the plant materials. CX, YZ and TX performed the research and analysed the data. CX and HZ wrote the paper. MK, XZ and HD revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This study was funded by the Central forestry reform and development fund -Representative National Key Protected Wild Plants Rescue and Protection in the Upper Reaches of the Yangtze River (zlg2021-cq20211210), and the Central forestry reform and development fund -Chengdu-Chongqing economic circle Joint Protection National Key Protected Wild Plants Rescue and Protection (zlg2022-cq20220907).